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Abstract 

Despite recent interest in reconstructing neuronal networks, complete wiring diagrams on the level of 
individual synapses remain scarce and the insights into function they can provide remain unclear. Even 
for Caenorhabditis elegans, whose neuronal network is relatively small and stereotypical from animal to 
animal, published wiring diagrams are neither accurate nor complete and self-consistent. Using materials 
from White et al. and new electron micrographs we assemble whole, self-consistent gap junction and 
chemical synapse networks of hermaphrodite C. elegans. We propose a method to visualize the wiring 
diagram, which reflects network signal flow. We calculate statistical and topological properties of the 
network, such as degree distributions, synaptic multiplicities, and small-world properties, that help in 
understanding network signal propagation. We identify neurons that may play central roles in information 
processing and network motifs that could serve as functional modules of the network. We explore 
propagation of neuronal activity in response to sensory or artificial stimulation using linear systems 
theory and find several activity patterns that could serve as substrates of previously described behaviors. 
Finally, we analyze the interaction between the gap junction and the chemical synapse networks. Since 
several statistical properties of the C. elegans network, such as multiplicity and motif distributions are 
similar to those found in mammalian neocortex, they likely point to general principles of neuronal 
networks. The wiring diagram reported here can help in understanding the mechanistic basis of behavior 
by generating predictions about future experiments involving genetic perturbations, laser ablations, or 
monitoring propagation of neuronal activity in response to stimulation. 



Introduction 

Determining and examining base sequences in genomes (U, J2l has revolutionized molecular biology. 
Similarly, decoding and analyzing connectivity patterns among neurons in nervous systems, the aim of 
the emerging field of connectomics ||3)-||5), may make a major impact on neurobiology. Knowledge of 
connectivity wiring diagrams alone may not be sufficient to understand the function of nervous systems, 
but it is likely necessary. Yet because of the scarcity of reconstructed connectomes, their significance 
remains uncertain. 

The neuronal network of the nematode Caenorhabditis elegans is a logical model system for advancing 
the connectomics program. It is sufficiently small that it can be reconstructed and analyzed as a whole. The 
302 neurons in the hermaphrodite worm are identifiable and consistent across individuals [7]. Moreover 
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the connections between neurons, consisting of chemical synapses and gap junctions, are stereotypical 
from animal to animal with more than 75% reproducibility lITTl- lPTOll . 

Despite a century of investigation ifTTTl . lfl2ll . knowledge of nematode neuronal networks is incomplete. 
The basic structure of the C. elegans nervous system had been reconstructed using electron micro- 
graphs Q, but a major gap in the connectivity of ventral cord neurons remained. Previous attempts 
to assemble the whole wiring diagram made unjustified assumptions that several reconstructed neurons 
were representative of others lfl3l . Much previous work analyzed the properties of the neuronal network 
(see e.g. fl4l - |[20l and references therein and thereto) based on these incomplete or inconsistent wiring 
diagrams Q, (T3). 

In this paper, we advance the experimental phase of the connectomics program Q, Ell by reporting 
a near-complete wiring diagram of C. elegans based on original data from White et al. Q but also 
including new serial section electron microscopy reconstructions and updates. Although this new wiring 
diagram has not been published before now, it has been freely shared with the community through the 
WormAtlas E2l and has also been used in studies such as E^l f*! 

We advance the theoretical phase of connectomics [24], |[25l , by characterizing signal propagation 
through the reported neuronal network and its relation to behavior. We compute for the first time, local 
properties that may play a computational purpose, such as the distribution of multiplicity and the number 
of terminals, as well as global network properties associated with the speed of signal propagation. Unlike 
the conventional "hypothesis-driven" mode of biological research, our work is primarily "hypothesis- 
generating" in the tradition of systems biology. 

Our results should help investigate the function of the C. elegans neuronal network in several ways. 
A full wiring diagram, especially when conveniently visualized using a method proposed here, helps 
in designing maximally informative optical ablation |[26l or genetic inactivation E71 experiments. Our 
eigenspectrum analysis characterizes the dynamics of neuronal activity in the network, which should help 
predict and interpret the results of experiments using sensory and artificial stimulation and imaging of 
neuronal activity. 

Organization of the RESULTS section reflects the duality of contribution and follows the tradition 
laid down by genome sequencing 12, (21 . We start by describing and visualizing the wiring diagram. 
Next, we analyze the non-directional gap junction network and the directional chemical synapse network 
separately. There are two primary reasons for separate analysis. First, understanding the parts before 
the whole provides didactic benefits. Second, separate consideration is valuable since we do not know 
the relative weight of gap junctions and chemical synapses and so any combination of the two involves 
additional assumptions. Finally, we analyze the combined network of gap junctions and chemical synapses. 

Results 

A. Reconstruction 

1) An Updated Wiring Diagram: The C. elegans nervous system contains 302 neurons and is divided 
into the pharyngeal nervous system containing 20 neurons and the somatic nervous system containing 
282 neurons. We updated the wiring diagram (see METHODS) of the larger somatic nervous system. 
Since neurons CANL/R and VC06 do not make synapses with other neurons, we restrict our attention 
to the remaining 279 somatic neurons. The wiring diagram consists of 6393 chemical synapses, 890 gap 
junctions, and 1410 neuromuscular junctions. 

The new version of the wiring diagram incorporates original data from White et al. [7], Hall and 
Russell ifTOl , updates based upon later work (8], |[28l , as well as new reconstructions. Although neuronal 

'See METHODS section for details on freely obtaining the wiring diagram in electronic form. 
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circuitry in the head and tail was previously documented Q, |[T0l . the connection details for 58 motor 
neurons in the ventral cord of the worm were lacking. We compiled most of the missing data using 
original electron micrographs and handwritten notes from White and coworkers. The dorsal side of the 
worm around the midbody, however, was not previously documented. Using original thin worm sections 
of animal N2U prepared by White et al. Q, we generated new micrographs and reconstructed neurons 
with processes in this region. In total, over 3000 synaptic contacts, including chemical synapses, gap 
junctions, and neuromuscular junctions were either added or updated from the previous version of the 
C. elegans wiring diagram. 

From our compilation of wiring data, including new reconstructions of ventral cord motor neurons, we 
applied self-consistency criteria to isolate records with mismatched reciprocal records. The discrepancies 
were reconciled by checking against electron micrographs and the laboratory notebooks of White et 
al. Connections in the posterior region of the animal were also cross-referenced with reconstructions 
published by Hall and Russell [ 10]. Reconciliation involved 561 synapses for 108 neurons (49% chemical 
'sends,' 31% chemical 'receives,' and 20% electrical junctions). The current wiring diagram is considered 
self-consistent under the following criteria: 

1) A record of Neuron A sending a chemical synapse to Neuron B must be paired with a record of 
Neuron B receiving a chemical synapse from Neuron A. 

2) A record of gap junction between Neuron C and Neuron D must be paired with a separate record 
of gap junction between Neuron D and Neuron C. 

Although the updated wiring diagram represents a significant advance, it is only about 90% complete 
because of missing data and technical difficulties. Due to sparse sampling along lengths of the sublateral, 
canal-associated lateral, and midbody dorsal cords, about 5% of the total chemical synapses are missing, 
as concluded from antibody staining for synapses |29l . Many gap junctions are likely missing due to the 
difficulty in identifying them in electron micrographs using conventional fixation and imaging methods. 
Hopefully, application of high-pressure freezing techniques and electron tomography will help identify 
missing gap junctions [30]. Finally, it should be noted that this reconstruction combined partial imaging 
of three worms, with images for the posterior midbody being from the male N2Y. 

The basic qualitative properties of the updated C. elegans nervous system remain as reported previ- 
ously Q-O- Neurons are divided into 118 classes, based on morphology, dendritic specialization, and 
connectivity. Based on neuronal structural and functional properties, the classes can be divided into three 
categories: sensory neurons, interneurons, and motor neurons. Neurons known to respond to specific 
environmental conditions, either anatomically, by sensory ending location, or functionally, are classified 
as sensory neurons. They constitute about a third of neuron classes. Motor neurons are recognized by the 
presence of neuromuscular junctions. Interneurons are the remainder of the neuron classes and constitute 
about half of all classes. A few of the neurons could have dual classification, such as sensory/motor 
neurons. Some interneurons are much more important for developmental function than for function in 
the final neuronal network 11501 . 

The majority of sensory neuron and interneuron categories contain pairs of bilaterally symmetric 
neurons. Motor neurons along the body are organized in repeating groups whereas motor neurons in 
the head have four- or six-fold symmetry. A large fraction of neurons send long processes to the nerve 
ring in the circumpharyngeal region to make synapses with other neurons Q. 

The neurons in C. elegans are structurally simple: most neurons have one or two unbranched processes 
and form en passant synapses along them. Dendrites are recognized by being strictly "postsynaptic" or by 
containing a specialized sensory apparatus, such as amphid and phasmid sensory neurons. Interneurons 
lack clear dendritic specialization. It is interesting to note that a given worm neuron has connections 
with only about 15% of neurons with which it has physical contact Q, (8), a similar number to the 
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Fig. 1. Adjacency matrices for the gap junction network (blue circles) and the chemical synapse network (red points) with 
neurons grouped by category (sensory neurons, interneurons, motor neurons). Within each category, neurons are in anteroposterior 
order. Among chemical synapse connections, small points indicate less than 5 synaptic contacts, whereas large points indicate 
5 or more synaptic contacts. All gap junction connections are depicted in the same way, irrespective of number of gap junction 
contacts. 



connectivity fraction in other nervous systems |3ll , |32l . 

2) Wiring Diagram as Adjacency Matrices: In the remainder of the paper, we describe and analyze 
the connectivity of gap junction and chemical synapse networks of C. elegans neurons. Gap junctions are 
channels that provide electrical coupling between neurons, whereas chemical synapses use neurotrans- 
mitters to link neurons. The network of gap junctions and the network of chemical synapses are initially 
treated separately, with each represented by its own adjacency matrix, Figure Q] In an adjacency matrix 
A, the element in the zth row and jth column, aij, represents the total number of synaptic contacts from 
the ith neuron to the jth. If neurons are unconnected, the corresponding element of the adjacency matrix 
is zero. An adjacency matrix may be used due to self-consistency in the gathered data. 

Although gap junctions may have directionality, i.e. conduct current in only one direction, this has 
not been demonstrated in C. elegans. Even if directionality existed, such information cannot be extracted 
from electron micrographs. Thus we treat the gap junction network as an undirected network with a 
symmetric adjacency matrix. Weights in both a-ij and aj% represent the total number of gap junctions 
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between neurons i and j. 

Since chemical synapses possess clear directionality that can be extracted from electron micrographs, 
we represent the chemical network as a directed network with an asymmetric adjacency matrix. The 
elements of the adjacency matrix take nonnegative values, which reflect the number of synaptic contacts 
between corresponding neurons. Contacts are given equal weight, regardless of the apparent size of the 
synaptic apposition. We use nonnegative values for most of the paper because we cannot determine 
whether a synapse is excitatory, inhibitory, or modulatory from electron micrographs of C. elegans. For 
the linear systems analysis, we do however make a rough guess of the signs of synapses based on 
neurotransmitter gene expression data. 

Electron micrographs for C. elegans have a further limitation that causes some synaptic ambiguity. 
When a presynaptic terminal makes contact with two adjacent processes of different neurons (send_joint 
in Durbin's notation [El), it is not known which of these processes acts as a postsynaptic terminal; both 
might be involved. We count all polyadic synaptic connections. Polyadic connections are briefly revisited 
in the DISCUSSION. 

3) Visualization: Although statistical measures that we investigate later in this paper provide significant 
insights, they are no substitute to exploring detailed connectivity in the neuronal network. As the number 
of connections between neurons is large even for relatively simple networks, such analysis requires a 
convenient way to visualize the wiring diagram. Previously, various fragments of the wiring diagram 
were drawn to illustrate specific pathways (H, ll33l . ll34l . Here, we propose a method to visualize the 
whole wiring diagram in a way that reflects signal flow through the network as well as the closeness of 
neurons in the network, Figure [2] To this end, we use spectral network drawing techniques because they 
have certain optimality properties |[35l and aesthetic appeal. Next, we give an intuitive description of our 
visualization method; mathematical details can be found in Appendix [A] 

The vertical axis in Figure Oa), represents the position of neurons in the signal flow hierarchy 1136*1 , 
lf37l of the chemical synapse network with sensory neurons at the top and motor neurons at the bottom, 
with interneurons in between. We want the vertical coordinate of pre- and post-synaptic neurons to differ 
by one, however due to "frustration" this is not always possible. Frustration happens when distances 
measured along network connections cannot be made to correspond to the hierarchy distances: there are 
two different hierarchical paths that require a particular neuron to appear in two different places. We look 
for the layout that has smallest deviation from this condition and find a closed form solution lf36l . [38]. 
The number of synapses from sensory to motor neurons — the signal flow depth of the network — can be 
read off the vertical coordinate. Depending on the specific neurons considered, the depth is typically 2-3 

ED. 

Neuronal position on the horizontal plane, Figure Ob), represents the connectivity closeness of neurons 
in the combined chemical and electrical synapse network. Neuronal coordinates are given by the second 
and third eigenmodes of the symmetrized network's graph Laplacian (see below). In this representation, 
pairs of synaptically coupled neurons with larger number of connections in parallel tend to be positioned 
closer in space. 

Thus, Figure [2] represents not the physical placement of neurons in the worm but signal flow and 
closeness in the network. Such visualization reveals that motorneurons and some interneurons segregate 
into two lobes along the first horizontal axis: the right lobe contains motorneurons in the ventral cord and 
the left lobe consists of neck/tail neurons. The bi-lobe structure suggests partial autonomy of motorneurons 
in the ventral cord and neck/tail. Interneurons that could coordinate the function of the two lobes can be 
easily identified by their central location. 
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Fig. 2. The C. elegans wiring diagram is a network of identifiable, labeled neurons connected by chemical and electrical 
synapses. Red, sensory neurons; blue, interneurons; green, motorneurons. (a). Signal flow view shows neurons arranged so that 
the direction of signal flow is mostly downward, (b). Affinity view shows structure in the horizontal plane reflecting weighted 
non-directional adjacency of neurons in the network. 
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B. Gap Junction Network 

For quantitative characterization, we first consider the gap junction network. 

1 ) Basic Structure and Connectivity: The gap junction network that we analyze consists of 279 neurons 
and 514 gap junction connections, consisting of one or more junctions. The network is not fully connected, 
but is divided into a giant component containing 248 neurons, two smaller components of 2 and 3 neurons, 
and 26 isolated neurons with no gap junctions (Table ISTT ). The giant component has 511 connections. 
An Erdos-Renyi random networlo with 279 neurons and connection probability 0.0133 (thus with 514 
expected connections) would be expected to have 271 neurons in the giant component. The true gap 
junction giant component is much smaller; the probability of finding such a small giant component in a 
random network is on the order of 10 -14 (see METHODS). A better comparison, however, can be made to 
random networks with degree distributions that match the degree distribution of the gap junction network 
j39l . Here, the degree of a neuron is the number of neurons with which it makes a gap junction. The giant 
component in a degree-matched random network would be expected to be 251 neurons (see METHODS), 
about the same size as the measured giant component. Using connectivity data from |fT3l , Majewska and 
Yuste had previously pointed out that most neurons in C. elegans belong to the giant component fiOl . 
Our results agree roughly with [40], although our dataset excludes non-neuronal cells and places certain 
neurons in different connected components. 

The adjacency matrix of the network, A, is depicted in Figure Q] (the number of gap junctions in a 
connection is not depicted). The matrix is symmetric since the network is undirected. We may explore 
the utility of representing the wiring diagram as a three-layer network by grouping neurons by category 
(sensory neurons, interneurons, motor neurons). As shown in Tables |S2"R and IS2B, each category has 
many recurrent connections; with the exception of connections between sensory and motor neurons, there 
are also many connections between categories. In particular, Table IS2B indicates that motor neurons send 
to interneurons roughly the same number of connections as recurrently sent back to motor neurons. These 
observations suggest that on the level of gap junctions, the value of a three-layer network abstraction is 
questionable. 

2) Distributions of Degree, Multiplicity and the Number of Terminals: In this section, we analyze sta- 
tistical properties of individual neurons and synaptic connections. To characterize the ability of individual 
neurons to propagate or collect signals, we compute the degree di of neuron i, which is the number of 
neurons that are coupled to i by at least one gap junction. The mean degree is 3.68, however this value 
is not representative as the degree varies in a wide range, from to 40. Thus, it is important to look 
at the degree distribution, which has been used to characterize and classify other networks previously 

ED-E3- 

To visualize the discrete degree distribution, p(d), we use the survival function: 

oo 

P(d) = 5>(*0. (D 

k=d 



which is the complement of the cumulative distribution function, Figure |3(a)[ The advantages of looking 
at the survival function rather than the degree distribution directly are that histogram binning is not 
required and that noise in the tail is reduced B31 . The survival function is also later applied to visualize 
other statistics. Various commonly encountered distributions and their corresponding survival functions 
are given in Appendix [B] 

2 Construction of an (unweighted) Erdos-Renyi random network requires a single parameter, the probability of a connection 
between two neurons. 



7 



We perform a fitting procedure for the tail of the gap junction degree distribution iPBI (see METHODS). 
We find that the tail (d > 4) can be fit by the power law with exponent 7 = 3.14, Figure [3(a)] but not 
by the exponential decay (p-value < 0.1). This result is consistent with the view that the gap junction 
network is scale-free [42 J. 

To characterize the direct impact that one neuron can have on another, we quantify the strength of 
connections by the multiplicity, rriij, between neurons i and j, which is the number of synaptic contacts 
(here gap junctions) connecting i to j. The degree treats synaptic connections as binary, whereas the 
multiplicity quantifies the number of contacts. The multiplicity distribution for the gap junction network 
is shown in Figure |3(b)| We find that the multiplicity distribution for m > 2 obeys a power law with 
exponent 7 = 2.76. Although the exponential decay fit to the tail passes the p- value test, the log-likelihood 
is significantly lower than for the power law. 

Finally, the number of terminals that lie on a given neuron is the sum of the multiplicities of all gap 
junction connections. The tail of the distribution of the number of synaptic terminals, Figure 3(c) is 
adequately fit by a power law with exponent 7 = 2.53. 




Fig. 3. Survival functions for the distributions of degree, multiplicity, and number of synaptic terminals in the gap junction 
network. Neurons or connections with exceptionally high statistics are labeled. The tails of the distributions can be fit by a 
power law with the exponent 3.14 for the degrees [(a)! 2.76 for the multiplicity distribution |(b)[ 2.53 for the number of synaptic 
terminals |(c)| The exponents for the power law fits of the corresponding survival functions are obtained by subtracting one. 



Identifying neurons that play a central or special role in the transmission or processing of information 
may also prove useful |46l- |[50l . To rank neurons according to their roles, we introduce several centrality 
indices. Perhaps the simplest centrality index is degree centrality Cd(i)- Degree centrality is simply the 
degree of a neuron, Cd(i) = rfj, and is motivated by the idea that a neuron with connections to many 
other neurons has a more important or more central role in the network than a neuron connected to only 
a few other neurons. Neurons that have unusually high degree centrality include AVAL/R and AVBL/R. 
The same neurons lie in the tail of the distribution of the number of synaptic terminals, Figure |3(c)[ 
suggesting strong coupling to the network. These neuron pairs are command interneurons responsible for 
coordinating backward and forward locomotion, respectively [22], |[34l , I5T1 . The high degree centralities 
of RIBL/R suggest a similarly central function for those neurons, though they each only have 19 gap 
junction terminals, in the middle of the distribution of number of terminals, suggesting weaker coupling 
to the network. 

3) Small World Properties: Having described statistical properties of individual neurons and connec- 
tions, such as the degree and multiplicity distributions, we now investigate properties that may describe 
the efficiency of signal transmission across the gap junction network. Traditionally [14], this analysis 
does not consider multiplicity of gap junctions but treats them as binary. We analyze signal propagation 
when including multiplicities in the next subsection. 
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The geodesic distance, dij, between two neurons in the network is the length of the shortest network 
path between them. The network path is measured by the number of connections that are crossed rather 
than by physical distance. The average geodesic distance over all pairs of neurons is the characteristic 
path length lfl4l : 

L = Niihr) £ (2) 

where N is the number of neurons. This global measure describes how readily or rapidly a signal can 
travel from one neuron to another since it is simply the average distance between all neurons. Clearly, 
the measure L requires the network to be connected (otherwise L diverges), so we restrict attention to 
the giant component. 

A signal originating in one neuron in the giant component must cross L = 4.52 gap junction connections 
on average to reach another neuron of the giant component. For an Erdos-Renyi random network with 248 
neurons and 511 connections, the characteristic path length is approximately log(248)/log(511/248) = 
7.63 lfl4l . When the actual degree distribution of the gap junction network is taken into account, a random 
network from that ensemble would be expected to have characteristic path length 3.05 (see METHODS). 



The distribution of geodesic distances dij in the giant component is shown in Figure S |l(a) 

A second measure for signal propagation is the clustering coefficient C, which measures the density 
of connections among an average neuron's neighbors. It is defined as lfT4l : 

where E{Mi) is the number of connections between neighbors of i, ki is the number of neighbors of 
i, and Ci measures the density of connections in the neighborhood of neuron i (we set d = 1 when 
ki = 1). We find the clustering coefficient C = 0.21. The clustering coefficient for an Erdos-Renyi 
random network with 248 neurons and 511 connections is approximately 511/(248 2 ) = 0.0083 lfl4~l . 
For a degree-matched random network, we computed the clustering coefficient C = 0.05. Thus, the 
giant component of the gap junction network is strongly clustered relative to random networks, both 
Erdos-Renyi and degree-matched. 

Small world networks are more clustered than Erdos-Renyi random networks and yet have smaller 
average distances [52]. Thus, the giant component of the gap junction network may be classified as a 
small world network. Table [S3] shows a comparison of the gap junction network of C. elegans with other 
networks that have been characterized as small world networks. 

Next we consider how quickly individual neurons reach all other neurons in the network. The normal- 
ized closeness of a neuron i is the average geodesic distance dij across all neurons j that are reachable 
from i |47l : 

^avg(i) = ^— j- di r ^ 

The normalized closeness centrality, which takes higher values for more central neurons, is defined as 
the inverse, c c (i) = l/cZ avg (i). 

Restricting to the giant component of the gap junction network, the six most central neurons are AVAL, 
AVBR, RIGL, AVBL, RIBL, and AVKL. In addition to command interneuron classes AVA and AVB, 
these include RIBL and RIGL, both ring interneurons, and AVKL, an interneuron in the ventral ganglion 
of the head. The set of neurons that are closeness central mostly overlaps with the set of neurons that 
are degree central. 
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The Spearman rank correlation coefficient ll53l between degree centrality Cd(i) and closeness centrality 
c c (i) for the entire giant component, however, is only 0.036. Since correlation between the two centrality 
measures does not extend to peripheral neurons, ordering of importance is different. 

4 ) Spectral Properties: Global network properties discussed in the previous section characterize signal 
transmission while ignoring connection weights. As weights affect the effectiveness of signal transmission 
and vary among connections, we now analyze the weighted network by using linear systems theory. 
Although neuronal dynamics can be nonlinear, spectral properties nevertheless provide important insights 
into function. For example, the initial success of the Google search engine is largely attributed to linear 
spectral analysis of the World Wide Web 11541 . 

We characterize the dynamics of the gap junction network by the following system of linear differential 
equations, which follow from charge conservation ||55l , [56]: 



where Vi is the membrane potential of neuron i, C{ is the membrane capacitance of neuron i, gij is the 
conductance of gap junctions between neurons i and j, and g™ is the membrane conductance of neuron i. 
Assuming that each neuron has the same capacitance C and each gap junction has the same conductance 
g, i.e. g^ = gAij, we can rewrite this equation in terms of the time constant r = C/g as: 



3 

Assuming that gap junction conductance is greater than the membrane conductance, we temporarily 
neglect the last term and rewrite this equation in matrix form: 

r% = ~LV, (7) 

where L is the Laplacian matrix of the weighted network, L = D — A, D contains the number of neuron 
gap junctions on the diagonal and zeros elsewhere, and V is a column vector of the membrane potentials. 
A different plausible differential equation model is discussed in Appendix O 

This system of coupled linear differential equations can be solved by performing a coordinate trans- 
formation to the Laplacian eigenmodes. Since the Laplacian eigenmodes are decoupled and evolve 
independently in time, performing an eigendecomposition of initial conditions leads to a full description of 
the system dynamics. We show the survival function of the eigenspectrum of the Laplacian in Figure [4(a) 

What insight can be gained from inspection of the Laplacian eigenmodes? The gap junction network 
is equivalent to a network of resistors, where each gap junction acts as a resistor. The eigenmodes give 
intuition about experiments where a charge is distributed among neurons of the network and the spreading 
charge among the neurons is monitored in time. If the charge is distributed among neurons according 
to an eigenmode, the relative shape of the distribution does not change in time. The charge magnitude 
decays with a time constant specified by the eigenvalue. The smallest eigenvalue of the Laplacian is 
always zero, corresponding to the infinite relaxation time. In the corresponding eigenmode each neuron 
is charged equally. 

If the charge is distributed according to eigenmodes corresponding to small eigenvalues, the decay is 
rather slow. Thus, these eigenmodes correspond to long-lived excitation. The existence of slowly decaying 
modes often indicates that the network contains weakly coupled subnetworks, in which neurons are 
strongly coupled among themselves. The corresponding charge distribution usually has negative values on 
one subnetwork and positive values on the other subnetwork. Because of the relatively slow equilibration 
of charge between the subnetworks, such eigenmode decays slowly. 




(5) 
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(6) 
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For example, one might speculate that the eigenmode associated with A3 (Figure 4(c) ) on the 'black' 
side reflects a coupling of chemosensory neurons in the tail (PHBL/R) along with interneurons (AVHL/R, 
AVFL/R) and motor neurons (VC01-05) involved in egg-laying behavior. At the level of gap junctions, 
these neurons are weakly coupled with chemosensory neurons in the head (ADFR, ASIL/R, AWAL/R) 
and related interneurons (AIAL/R) on the 'red' side. 

Another interesting example is the eigenmode associated with A 13 (Figure |4(d)| ). Neurons on the 'red' 
side overlap significantly with those identified previously in a hub-and-spoke circuit mediating pheromone 
attraction, oxygen sensing, and social behavior 1571 . Such overlap is consistent with the view |57] that 
this network of neurons solves a consensus problem [58]. 

The above two examples demonstrate that spectral analysis can uncover circuits that have been described 
using experimental studies. The probability of a known functional circuit to appear in an eigenmode by 
chance is small (see METHODS). It would be interesting to see whether other eigenmodes have a biological 
interpretation and therefore generate predictions for future experiments. 

To prioritize further analysis of eigenmodes for biological significance, it may be advantageous to 
focus on the slow and sparse modes, where few neurons exhibit significant activity. We can quantify 
sparseness of normalized eigenmodes by the sum of absolute values (rectilinear norm) of the eigenmode 
components; sparser eigenmodes have smaller rectilinear norms |59l . Figure |4(b)| is a scatterplot of 
eigenmodes showing both their decay constant and their rectilinear norms. 

The full set of eigenmodes of the connected component is shown in Figure [S2] The eigenmodes 
corresponding to large eigenvalues decay fast, suggesting that corresponding neurons have the same 
membrane potential on relevant time scales and act effectively as a single unit. Many such eigenmodes 
peak (with opposite signs) for left-right neuronal pairs (Figure |S3l ), often known to be functionally 
identical, which therefore act as a single unit. 

What is the absolute value of decay constants for various eigenmodes? Current knowledge of electrical 
parameters for C. elegans neurons allows us to estimate the decay times only approximately. Assuming 
neuron capacitance of 2pF [60] and gap junction conductance of 200pS, we find a time constant r = 10ms. 
This implies that the slowest non-trivial mode corresponding to the second lowest eigenvalue, A2 = 0.12 
has decay time of about 83ms, Figure 4(a)| This eigenvalue, A2, is known as the algebraic connectivity 
of a network [61] and is discussed further in Appendix ICl 

What is the effect of the dropped term corresponding to the membrane current in ©? As this term 
would correspond to adding a scaled identity matrix to the Laplacian, the spectrum should uniformly 
shift to higher values by the corresponding amount. Thus, even the eigenmode corresponding to the zero 
eigenvalue would now have a finite decay time. Assuming the membrane conductance of about lOOpS 
[60], we find 20ms decay time. This leads to a 0.5 increase in the values of A. Now, the slowest non-trivial 
mode corresponds to a decay time of about 16ms. 

In addition to highlighting groups of neurons that could be functionally related, spectral analysis allows 
us to predict, under linear approximation, the outcome of experiments that study the spread of an arbitrarily 
generated excitation in the neuronal network. Such excitation can be generated in sensory neurons by 
presenting a sensory stimulus ll62l or in any neuron by expressing channelrhodopsin in that cell |]26l . ll63l . 
|[64l . The spread of activity can be monitored electrophysiologically or using calcium-sensitive indicators. 

To predict the spread of activity, we may decompose the excitation pattern into the eigenmodes and, 
by taking advantage of eigenmode independence, express temporal evolution as a superposition of the 
independently decaying eigenmodes. The initial redistribution of charge would correspond to the fast 
eigenmodes, whereas the long-term evolution of charge distribution would be described by the slow 
eigenmodes. Appendix [D] further discusses eigendecomposition and the interpretation of eigenmodes. 
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(c) (d) 

Fig. 4. Linear systems analysis of the giant component of the gap junction network. |(a)| Survival function of the eigenvalue 
spectrum (blue). The algebraic connectivity, A2, is 0.12 and the spectral radius, A248, is 118. A time scale associated with 
the decay constant is also given. |(b)| Scatterplot showing the rectilinear norm and decay constant of the eigenmodes of the 
Laplacian. The fastest modes from Figure [S3] are marked in red. The sparsest and slowest modes, most amenable to biological 
analysis, are located in the lower-left corner of the diagram. |(c)[ Eigenmode of Laplacian corresponding to A3 (marked green 
in panel [(b)} . |(d)| Eigenmode of Laplacian corresponding to A13 (marked cyan in panel [(b)| >. 



5) Motifs: Several of the quantitative properties computed thus far measure global network structure 
and may determine aspects of system operation. Now we examine the network locally and analyze the 
frequency of various connectivity subnetworks among small groups of neurons. Overrepresentation in the 
subnetwork distribution often displays building blocks of the network such as computational units [17], 
||65l . Since the gap junction network is undirected, there are four kinds of subnetworks that can appear 
over three neurons; this distribution is shown in Figure 5(a)[ As a null-hypothesis we use random network 
ensembles that preserve the degree distribution. We find that fully connected triplets are overrepresented. 

Four neurons can be wired into 11 kinds of subnetworks; this distribution is shown in Figure [5(b)] 
In the case of quadruplets, the null-hypothesis preserves the degree for each neuron and the number 
of triangles. A numerical rewiring procedure is used to generate samples from these random network 
ensembles ||39l , 1166*1 . since no analytical expression for expected subnetwork counts is extant 116711 . We 
find that a "fan" (motif #7) and a "diamond" (motif #10) are overrepresented. 
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Note that neurons participating in motifs also make connections with neurons outside of the motif, 
which are traditionally not drawn in putative functional circuits JSJ, Il62l . Such putative functional circuit 
diagrams may even omit connections within the motif (SJ, (62), which we do not allow. 
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Fig. 5. Subnetwork distributions for the gap junction network. Overrepresented subnetworks are boxed, with the p- value from 
the step-down min-P-based algorithm for multiple-hypothesis correction |66|, |68| (n = 1000) shown inside, [(a)] The ratio of 
the 3-subnetwork distribution and for the mean of a degree-preserving ensemble of random networks (n — 1000). The counts 
for the particular random networks that appeared in the ensemble are also shown. |(b)| The ratio of the 4-subnetwork distribution 
and for the mean of a degree and triangle-preserving ensemble of random networks (n = 1000). The counts for the particular 
random networks that appeared in the ensemble are also shown. 



C. Chemical Synapse Network 

Now we consider the chemical synapse network. Recall that due to structural differences between 
presynaptic and postsynaptic ends of a chemical synapse, electron micrographs can be used to determine 
the directionality of connections. Hence the adjacency matrix is not symmetric as it was for the gap 
junction network. 

1) Basic Structure and Connectivity: The network that we analyze consists of 279 neurons and 2194 
directed connections implemented by one or more chemical synapses. The adjacency matrix of the 
network shown in Figured] is suggestive of a three-layer architecture. Table |S4l shows the distribution of 
connections between categories in the three-layer architecture. Each chemical subnetwork is characterized 
by a high number of recurrent connections, just as for the gap junction. However, the majority of 
connections with other subnetworks is consistent with feedforward information processing (sensory to 
interneuron and interneuron to motorneurons). Therefore, a three-layer network abstraction may be more 
valuable for chemical synapses than for gap junctions. 

There are two different definitions of connectivity for directed networks. A weakly connected compo- 
nent is a maximal group of neurons which are mutually reachable by possibly violating the connection 
directions, whereas a strongly connected component is a maximal group of neurons that are mutually 
reachable without violating the connection directions. The whole chemical synapse network is weakly 
connected and can be divided into a giant strongly connected component with 237 neurons, a smaller 
strongly connected component of 2 neurons, and 40 neurons that are not strongly connected (Table |S5l ). 
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The random directed network corresponding to the chemical network is fully weakly connected, even 
when the degree distribution is taken into account (see METHODS). A strongly connected giant component 
as small as in the chemical network is not likely in a random network (see f69l ). Thus, the chemical 
network is more segregated than would be expected for a random network. 

2) Distributions of Degree, Multiplicity and the Number of Terminals: Since chemical synapses form a 
directed network, neuron connectivity is characterized by in-degrees (the number of incoming connections) 
and out-degrees (the number of outgoing connections) rather than simply degrees. The joint distribution 
of in-degrees and out-degrees is shown in Figure 6(a) As can be seen by the distribution clustering 
around the diagonal line, the in-degrees and out-degrees are correlated; the correlation coefficient is 0.52, 
very similar to the correlation coefficient of email networks, 0.53 |[70l . 

The survival functions associated with the marginal distributions of in-degrees and out-degrees are 
shown in Figures [6(b)| and [6(c)| respectively The mean number of incoming and outgoing connections is 
7.86 each. We attempt to fit these distributions. The tails of the two distributions can be satisfactorily fit 
by power laws with exponents 3.17 and 4.22 respectively. Exponential fit is ruled out (p-value < 0.1) for 
the in-degree but not for the out-degree distribution. In the latter case, the log-likelihood is insignificantly 
lower for the exponential decay than for the power law. 

Multiplicity of connection, m^, is the number of synapses in parallel from neuron i to neuron j. The 
corresponding survival function (including unconnected pairs) is shown in Figure |6(d)| The mean number 
of synapses per connection (excluding unconnected pairs) is 2.91. The tail of the distribution can be fitted 
by an exponential, but not by a power law (p-value < 0.1). In addition, the whole distribution (m > 1) 
can be fit by a stretched exponential (or Weibull) distribution, p(m) ~ (m//3) 7 ~ 1 e _ ^ m// ^^ with the scale 
parameter f3 = 0.36 and the shape parameter 7 = 0.47. A stretched exponential applied to the whole 
distribution has the same number of fitting parameters as an exponential decay fit to the tail starting with 
an adjustable m. Log-likelihood comparison of the tail exponential and the stretched exponential favors 
the latter. 

As for the gap junction network, we can also study the distribution of number of synaptic terminals on 
a neuron. This involves adding the multiplicities of the connections, rather than just counting the number 
of pre- or post-synaptic partners. The joint histogram (not shown) exhibits similar correlation as for the 
degree distribution, with correlation coefficient 0.42. 

Figures 6(e) and 6(f) show the marginal survival functions for the number of post-synaptic terminals 
(in-number) and the number of pre-synaptic terminals (out-number). The mean number of pre- and 
post-synaptic terminals is 22.9 each. We were unable to find a satisfactory simple fit to the in-number 
distribution, Figure [6(e)| The tail of the out-number distribution could be fit by a power law with exponent 
4.05, but not by an exponential, Figure 6(f)| 

As for the gap junction network, we can identify central neurons (cf. |49l , fTTI ) for the chemical 
network. The degree centrality in a directed network may be defined with respect to the in-degree or the 
out-degree. Interestingly, neuron AVAL has the best in-degree, whereas AVAR has the best out-degree 



and AVAR has the best out-degree and AVAL has the second best out-degree, Figure 6(a) 



3) Small World Properties: In the strongly connected component, we can define the directed geodesic 
distance as the shortest path between two neurons that respects the direction of the connections. This 
distribution is shown in Figure S |l(b)| The directed characteristic path length, L, is the average directed 
geodesic distance over all pairs of neurons in the strongly connected component and is computed to be 
3.48 steps. For a random network degree-matched to the chemical network, one would expect L = 2.91. 
The similarity of the geodesic distances suggests that signals diffuse as quickly as in a random network. 

Although there are several definitions of clustering for directed graphs in the literature 11721 . we use the 
clustering of the out-connected neighbors since it captures signal flow emanating from a given neuron. 
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Fig. 6. Degree distribution |(a)| and survival functions for the distributions of in-/out-degree, multiplicity, and in-/out-number of 
synaptic terminals in the chemical synapse network [(b)(|(f)| Neurons or connections with unusually high statistics are labeled. 
The tails of the distributions can be fit by a power law with exponents 3.17 for in-degree |(b)[ 4.22 for out-degree |(c)| and 4.05 
for out-number |(f) | The exponents for the survival function fits can be obtained by subtracting one. The survival function of the 
multiplicity distribution for m > 1 can be fit by a stretched exponential of the form e - ^ 71 ^ where j3 = 0.36 and 7 = 0.47 
|(d)| No satisfactory fit was found for the distribution of in-numbers |(e)[ 
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where E{Mi) is the number of connections between out-neighbors of neuron i, k{ is the number of 
out-neighbors of i, and Cj measures the density of connections in the neighborhood of neuron i. For the 
chemical network, the clustering coefficient is 0.22. Using the Watts-Strogatz approximations to L and 
C, the clustering coefficient for a random network is C r « i exp( ln ^ ); so for iV = 279 and L = 3.48, 
a random network would have C ~ 0.018. For a degree-matched random network we computed the 
clustering coefficient C = 0.079. Since the clustering coefficient for the chemical network is much more 
than a similar random directed network, it may be considered a small-world network, cf . Table [S3] 

For directed networks, measures of in-closeness and out-closeness may be defined using the average 
directed geodesic distance. In particular, the normalized in-closeness is the average geodesic distance 
from all other neurons to a given neuron: 

dfavgW = Yl d ^ (9) 

and the out-closeness is the average geodesic distance from a given neuron to all other neurons: 

where N is the number of neurons. Normalized centralities are the inverses: Cj C (i) = l/di avg (i) and 
c c(i) = l/^oavg(i)- The motivation behind these indices is similar to that in the gap junction case. 
In-closeness central neurons can be easily reached from all other neurons in the network. Out-closeness 
central neurons can easily reach all other neurons in the network. Normalized in-closeness centrality Ci C (i) 
and normalized out-closeness centrality c oc (i) are weakly anti-correlated, with correlation coefficient 
-0.12. 

For the giant component of the chemical network, the most in-closeness central neurons include AVAL, 
AVAR, AVBR, AVEL, AVER, and AVBL. All are command interneurons involved in the locomotory 
circuit; these neurons are also central in the gap junction network. The in-closeness centrality of command 
interneurons may indicate that in the C. elegans nervous system, signals can propagate efficiently from 
various sources towards these neurons and that they are in a good position to integrate it. 

The most out-closeness central neurons include DVA, ADEL, ADER, PVPR, AVJL, HSNR, PVCL, 
and BDUR. Only PVCL is a command interneuron involved in locomotion. The neuron DVA is an 
interneuron that performs mechanosensory integration; ADEL/R are sensory dopaminergic neurons in 
the head; and the other central neurons are interneurons in several parts of the worm. The out-closeness 
centrality of these neurons may indicate that signals can propagate efficiently from these neurons to the 
rest of the network and that they are in a good position for broadcast. 

4) Spectral Properties: Although chemical synapses are likely to introduce more nonlinearities than 
gap junctions, linear systems analysis can provide interesting insights, especially in the absence of 
other tools. Such an approach has additional merit in C. elegans, where neurons do not fire classical 
action potentials ll60l and have chemical synapses that likely release neurotransmitters tonically (56). To 
justify such analysis, a system of linear equations may be derived by approximating sigmoidal synaptic 
transmission functions with linear dependencies. This can be done by expanding synaptic transmission 
functions into a Taylor series around an equilibrium point ll56l . 

A major source of uncertainty in linear systems analysis of the chemical network is the unknown 
sign of connections, i.e. excitatory or inhibitory, due to the difficulty in performing electrophysiology 
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experiments. We use a rough approximation that GABAergic synapses are inhibitory, whereas glutamergic 
and cholinergic synapses are excitatory [73 ], but see Il6*2l . Thus inhibitory neurons are identified by looking 
at GAB A expression f74ll § 

Similarly to the gap junction network, we write the system of linear differential equations for the 
chemical synapse network [55], 1156*1 : 

Ci^ = J2 V Mi-9rVi, (11) 

j 

where V{ is the membrane potential of neuron i measured relative to the equilibrium, C{ is the membrane 
capacitance of neuron i, gji is the conductance in neuron i contributed by a chemical synapse in response 
to voltage Vj measured relative to the equilibrium and gf 1 is the membrane conductance of neuron i. 
Assuming that each neuron has the same capacitance C and each chemical synapse contact has the same 
conductance g, i.e. g±j = gA^j, we can rewrite this equation in terms of the time constant r = C/g as: 

= X>i A * " f V «- (12) 

j 

To avoid redundancy we defer analyzing this system of differential equations to the next section, where 
we consider the combined system including both gap junctions and chemical synapses. 

5) Motifs: We also find subnetwork distributions for the chemical synapse network. Since the network 
is directed, there are many more possible subnetworks. In particular there are 3 possible subnetworks on 
two neurons and 16 possible subnetworks on three neurons. We identify overrepresented subnetworks by 
comparing to random networks, generated with a rewiring procedure |f3~9l , ll66l . Such random network 
ensembles preserve in-degree and out-degree in the case of doublets and, additionally, the numbers of 
bidirectional and unidirectional connections for each neuron in the case of triplets. 

Figures 7(a) and |7(b)| show the subnetwork distributions on two and three neurons, respectively. We 



find that the C. elegans network contains similar overrepresented subnetworks as found by analyzing 
incomplete data |fT71 , [66]. For example, there is greater reciprocity in the chemical network than would 
be expected in a random network. Similarly, triplets with connections (of any direction) between each 
pair of neurons (seven rightmost triplets in Figure |7(b)[ ) collectively occur with much greater frequency 
than would be expected for a random network. 

Overrepresentation of reciprocal [8, Ch. 7] and triangle [7] motifs were previously noted. Such over- 
representation would arise naturally if proximity was a limiting factor for connectivity, however there is 
no evidence for this limitation. Rather we believe motifs have a functional role. 



D. Full Network 

Having considered the gap junction network and the chemical synapse network separately, we also 
examine the two networks collectively. To study the two networks, one may either look at a single network 
that takes the union of the connections of the two networks or one may look at the interaction between 
the two networks. 

1 ) Single Combined Network: First we look at a combined network, which is produced by simply 
adding the adjacency matrices of the gap junction and chemical networks together, while ignoring con- 
nection weights. Thus we implicitly treat gap junction connections as double-sided directed connections. 
This new network consists of 279 neurons and 2990 directed connections. It has one large strongly 
connected component of 274 neurons and 5 strongly isolated neurons. The five isolated neurons are 

3 The 26 GABAergic neurons are DVB, AVL, RIS, DD01-DD06, VD01-VD13, and the four RME neurons. 
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Fig. 7. Subnetwork distributions for the chemical synapse network. Overrepresented subnetworks are boxed, with the p-value 
from the step-down min-P-based algorithm for multiple-hypothesis correction |66|, |68| (n = 1000) shown inside, [(a)] The 
ratio of the 2-subnetwork distribution and the mean of a random network ensemble (n = 1000). Realizations of the random 
network ensemble are also shown. |(b)| The ratio of the 3-subnetwork distribution and the mean of a random network ensemble 
(n — 1000). Realizations of the random network ensemble are also shown. 



IL2DL/R, PLNR, DD06, and PVDR; this set is simply the intersection of the isolated neurons in the 
gap junction and chemical networks and does not seem to have any commonalities among members. Of 
course, it follows that since the chemical network is a single weakly connected component, this combined 
network is also a single weakly connected component. 

Naturally, the combined network is more compact than the individual networks. The mean path length 
L = 2.87, the geodesic distance distribution (Figure £ 1(c) ) becomes narrower. For a random network 
degree-matched to the combined network, one would expect L = 2.62, not significantly different. The 
clustering coefficient for the combined network is C = 0.26. The clustering coefficient for a similar 
random network would have been C = 0.026 lfT4l . and for a degree-matched random network C = 0.10. 
Therefore, the combined network, just like the individual networks, may be classified as small world. 
Turning to closeness centrality, the most in-close central neurons are AVAL/R, AVBR/L, and AVEL/R, as 
would be expected from the individual networks. The most out-close central neurons are DVA, ADEL, 
AVAR, AVBL, and AVAL, which include the top out-close neurons for both individual networks. 

We can also calculate the degree distribution of this combined network. The correlation coefficient 
between the in-degree and out-degree is 0.71; it is not surprising that the coefficient is so large considering 
that the gap junctions introduce an in- and out-connection simultaneously. Similar to the chemical synapse 
network, the tails of both the in-degree and the out-degree survival functions (Figures S |4(a)| and S |4(b)| ) 
can be fit with power laws. The tail of the out-degree could also be fit by an exponential decay, albeit 
with lower likelihood. 

The neurons with the greatest degree centrality are AVAL and AVAR. As for the chemical synapse 
network, neuron AVAL has the best in-degree and AVAR has the second best in-degree, whereas AVAR 
has the best out-degree and AVAL has the second best out-degree (Figures S 4(a) and S |4(b)| ). The next 
two neurons are AVBL/R in both in-degree and out-degree senses. 

As for the chemical synapse network, the tail of the out-number distribution was fit by a power law and 
the tail of the in-number distribution could not be fit satisfactorily. The tail of the out-number distribution 
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could also be fit by an exponential, albeit with lower likelihood. The multiplicity can be fit satisfactorily 
by a stretched exponential. 

2) Spectral properties: In this section we apply linear systems analysis to the combined network of 
chemical synapses and gap junctions taking into account multiplicities of individual connections. Due to 
our ignorance about the relative conductance of a single gap junction and of a single chemical synapse, 
we assume that they are equal. By combining equations © and (fT2l we arrive at: 

t% = D-^Sf p) + ^ hem) ] - fvi- w 

where yl^ hem ^ is negative if neuron j is GABAergic and positive otherwise. 

We proceed to find a spectral decomposition for the combined network. To avoid trivial eigenmodes, 
we restrict our attention to the strongly connected component of the combined network containing 274 
neurons. As before, we ignore the term and only study the matrix = — L(s a P) + A T ( chem ). Since $ 
is not symmetric, eigenvalues and eigenmodes may be complex-valued, occurring in complex conjugate 
pairs. Eigenvalues are plotted in the complex plane in Figure |8(a) 



What is the meaning of complex eigenvalues? The imaginary part of an eigenvalue is the frequency 
at which the associated eigenmode oscillates. The real part of an eigenvalue determines the amplitude 
of the oscillation as it varies with time. Eigenmodes that have an eigenvalue with a negative real part 
decay with time, whereas eigenmodes that have an eigenvalue with a positive real part grow with time. 



one should keep in mind that the ignored ^-Vi term would shift the real part of the eigenvalues towards 



When examining the temporal evolution of the eigenmodes whose eigenvalues are shown in Figure |8(a) 

L111U L11CIL Lilt, lyiX"." — 

more negative values 

As for the gap junction network alone, we can look for eigenmodes that may have functional signif- 
icance. For example, the sixth eigenmode of the combined network, Figure |8(b)[ includes neurons that 
are involved in sinusoidal body movement. As before, one may focus on sparse and slow eigenmodes 
for ease of investigation. The distribution of rectilinear norm and real part of eigenvalues is shown in 
Figure |8(c)[ and twelve of the sparsest and slowest eigenmodes of the combined network are plotted in 
Figure [8(d)] 

Having the eigenspectrum of the combined network allows one to calculate the response of the network 
to various perturbations. By decomposing sensory stimulation among the eigenmodes and following 
the evolution of each eigenmode, one could predict the worm's response to the stimulation. A similar 
calculation could be done for artificial stimulation of the neuronal network, induced for example, using 
channelrhodopsin [261, I3J, US - 

3 ) Interaction Between Networks: We have measured the structural properties of the combined network 
formed by adding together the adjacency matrices of the gap junction and chemical synapse network, 
however it is unclear how they interact. The two networks could be independent, or their connections 
could overlap more or less often than by chance. 

To investigate how the two networks overlap, we look at local structure. Figure [9] shows the likelihood 
ratios of chemical synapse connections being absent, being unidirectional, and being bidirectional given 
the presence or absence of a gap junction between the same pair of neurons (see METHODS). As can 
be seen, chemical synapses are more likely to be absent when there is no gap junction than when there 
is one. Unidirectional, and especially bidirectional, chemical synapses are more likely when there is a 
gap junction between given neurons. In this sense, the two networks are correlated, however it should be 
noted that when there is a gap junction, about 60% of the time there is no chemical synapse in either 
direction either. 
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Fig. 8. Linear systems analysis for the strong giant component of the combined network . |(a)| Eigenvalues plotted in the complex 
plane, [(b)] The eigenmode associated with eigenvalue Ae (marked cyan in panel [(c)) . |(c)| Scatterplot showing the sparseness 
and decay constant of the eigenmodes. |(d)| Sparse and slow eigenmodes of the combined network (marked red in panel \(c)\ . 
The real parts of the eigenmodes corresponding to A222, A224, A225, A226, A227, A232, A267, A268, A270, A272, A273, and A274 are 
shown. The eigenmodes are labeled with neurons that take value above a fixed absolute value threshold. Neurons with negative 
values are in red, whereas neurons with positive values are in black. 



Durbin had found that chemical and gap junction networks are essentially independent when imposing 
physical adjacency restrictions [8, Ch. 7], but as noted above, there is no evidence that proximity is a 
limiting factor for connectivity. We believe there may be a functional role for correlation/anticorrelation 
of the joint presence of gap junction and chemical connections. 

Why might the presence of connections in two networks either be correlated or anticorrelated? One 
possibility is that correlated connections simultaneously perform different functions [75] whereas anti- 
correlation yields connections between distinct kinds of neuronal pairs ll76l -[78l. 

What are the different functions performed by chemical synapses and gap junctions that could lead 
to correlation? One possibility is that the two different functions are sign-inverting and non-inverting 
coupling. Gap junctions are non-inverting: higher potential in a neuron raises the potential in other gap- 
junction-coupled neurons. Chemical synapses, on the other hand, may be either excitatory (non-inverting) 
or inhibitory (inverting). When the likelihood computations are repeated considering only neuron pairs 
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TABLE I 

Degree Sequence Correlation Coefficients 





gap/in 


gap/out 


in/out 


email [ 70 1 


correlation coefficient p 


0.64 


0.44 


0.52 


0.53 


avg. rand. perm, p 


-0.00 ±0.06 


0.00 ± 0.06 


0.00 ± 0.06 





where the presynaptic neuron is known to be GABAergic |74l , there is not much change, see Figure [9] 
This suggests that the primary purpose of overlapping inhibitory chemical synapses is not to counter 
excitatory gap junctions. This result, however, is only suggestive since the neurotransmitters and their 
action on postsynaptic receptors in many neurons have not been determined. Some other reason, such as 
differing temporal properties or robustness from redundancy, is needed to explain correlation. 

Another measure of the interaction between the two networks is the correlation between the degree 
sequences. The correlation coefficient between the gap junction degree and the chemical network in-degree 
is greater than and the correlation coefficient between the gap junction degree and the chemical network 
out-degree is less than the correlation coefficient between the chemical network in-degree and out-degree, 
as shown in Table U where comparisons to correlation coefficients between randomly permuted degree 
sequences (see METHODS) are shown. Large correlation coefficients imply that neurons are ordered in 
similar ways according to degree centrality. 

The two networks seem to primarily reinforce each other with correlated structure rather than augment 
each other with anticorrelated connections. 
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Fig. 9. Likelihood ratio for chemical network doublets given the absence/presence of a gap junction between the two neurons. 



E. Robustness Analysis 

Although the reported wiring diagram corrects errors in previous work and is considered self-consistent, 
one might wonder how remaining ambiguities and errors in the wiring diagram might affect the quanti- 
tative results presented. For network properties that are defined locally, such as degree, multiplicity, and 
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subnetwork distributions, clearly small errors in the measured wiring diagram lead to small errors in the 
calculated properties. For global properties such as characteristic path length and eigenmodes, things are 
less clear. 

To study the robustness of global network properties to errors in the wiring diagram, we recalculate 
these properties in the wiring diagrams with simulated errors. We simulate errors by removing randomly 
chosen synaptic contacts with a certain probability and assigning them to a randomly chosen pair of 
neurons. Then, we calculate the global network properties on the ensemble of edited wiring diagrams. 
The variation of the properties in the ensemble gives us an idea of robustness. 

First, we explore the robustness of the small world properties and the giant component calculations. We 
edit wiring diagrams by moving each gap junction contact with 10% probability and chemical synapse 
contact with 5% probability. Tables |S6] and [S7] show the global properties for 1000 random networks 
obtained by editing the experimentally measured network. These tables suggest that our quantitative 
results are reasonably robust to ambiguities and errors in the wiring diagram. 

Properties for the neuronal network from prior work in [13 ] are also shown for comparison. The number 
of synaptic contacts that must be moved to achieve this network (editing distance) roughly corresponds 
to that with 25.6% probability. 

Second, we characterize robustness for the linear systems analysis. Because of greater sensitivity of the 
eigenvalues to errors, we edit wiring diagrams by moving each gap junction contact with 1% probability 
and a chemical synapse contact with 0.5% probability. The spectra for 100 randomly edited networks 
along with the spectrum for the measured network (Figure 8(a) I are shown in Figure [10] Although the 
locations of eigenvalues shift in the complex plane, many of them move less than the nearest neighbor 
distance and remain isolated. 

In addition to considering the effect of typical random edits, we can characterize the effect of worst- 
case errors on the eigenvalues using the e-pseudospectrum fl79l , which gives the eigenvalue loci A e for 
all perturbations by matrices of norm e (Figure [TOi For the gap junction, A € (L) is simply the set of disks 
of radius e around the eigenvalues, but for the chemical and combined networks, A e (A T ) and A e (<3?) are 
larger. In the worst case scenario, most eigenmodes become mixed up. 

Electron micrographs of chemical synapses have a further ambiguity when more than one postsynaptic 
partner receives input at a release site. We treated such polyadic (sendjoint) synapses no differently than 
other synapses, but one might alternatively determine multiplicity by counting such synapses at 50% 
strength. This alternate quantitation clearly does not change statistics that ignore multiplicity; the change 
in the spectrum is depicted in Figure [TOl 



Discussion 

We have presented a corrected and more comprehensive version of the neuronal wiring diagram of 
hermaphrodite C. elegans using materials from White et al. [7] and new electron micrographs. We also 
proposed a convenient way to visualize the neuronal wiring diagram. The corrected wiring diagram and 
its visualization should help in planning experiments, such as neuron ablation. 

Next, we performed several statistical analyses of the corrected wiring, which should help with inferring 
function from structure. 

By using several different centrality indices, we found central neurons, which may play a special role 
in information processing. In particular, command interneurons responsible for worm locomotion have 
high degree centrality in both chemical and gap junction networks. Interestingly, command interneurons 
are also central according to in-closeness, implying that they are in a good position to integrate signals. 
However, most command interneurons do not have highest out-closeness, meaning that other out-closeness 
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Fig. 10. The spectrum of the giant component of the combined network matrix <1? (red), e-disks around the spectrum (light 
blue), spectra of 100 randomly edited networks (blue), and the e-pseudospectrum (orange). The value e = 4 is used (the average 
spectral norm of the 100 editing matrices was 3.4 ±0.9). The spectrum of the giant component of the combined network matrix 
<!> under an alternate quantitation of send_Joint synapses is also shown (green). 



central neurons, such as DVA, ADEL/R, PVPR, etc., are in a good position to deliver signals to the rest 
of the network. 

Linear systems analysis yielded a principled methodology to hypothesize functional circuits and to 
predict the outcome of both sensory and artificial stimulation experiments. We have identified several 
modes that map onto previously identified behaviors. 

Networks with similar statistical structural properties may share functional properties thus providing 
insight into the function of the C. elegans nervous system. To enable comparison of the C. elegans 
network with other natural and technological networks [80], we computed several structural properties 
of the neuronal network. In particular, the gap junction network, the chemical synapse network, and the 
combined neuronal network may all be classified as small world networks because they simultaneously 
have small average path lengths and large clustering coefficients |[T4l . 

The tails of the degree and terminal number distributions for the gap, chemical and combined networks 
(with the exception of the in-numbers) follow a power law consistent with the network being scale-free 
in the sense of Barabasi and Albert [42]. The tails of some distributions can also be fit by an exponential 
decay, consistent with a previous report [fl~5l. However, we found that exponential fits for the tails have 
(sometimes insignificantly) lower log-likelihoods than power laws, making the exponential decay a less 
likely alternative. For whole distributions, neither distribution passes the p-value test; if one is forced to 
choose, the exponential decay may be a less poor alternative. 

Several statistical properties of the C. elegans network are similar to those of the mammalian cortex. 
In particular, the whole distribution of C. elegans chemical synapse multiplicity is well-fit by a stretched 
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exponential (or Weibull) distribution (Figure |6(d)| ). Taking multiplicity as a proxy of synaptic connec- 
tion strength, this is reminiscent of the synaptic strength distribution in mammalian cortex, which was 
measured electrophysiologically, ll32l . ll68l . The definition of stretched exponential distribution is slightly 
different ['32], but has the same tail behavior. The stretch factor is ~ 0.5, close to that in the cortical 
network. 

In addition, we found that motif frequencies in the chemical synapse network are similar to those in the 
mammalian cortex [ 68 ] . Both reciprocally connected neuron pairs and triplets with a connection between 
every pair of neurons (regardless of direction) are over-represented. The similarity of the connection 
strength and the motif distributions may reflect similar constraints in the two networks. Since proximity 
is unlikely to be the limiting factor, we suggest that these constraints may reflect functionality. We found 
that the chemical synapse and the gap junction networks are correlated, which may provide insight into 
their relative roles. 

To conclude the paper, let us note that our scientific development was not hypothesis-driven, but rather 
exploratory. Yet we hope that the reported statistics will help in formulating a theory that explains how 
function arises from structure. 

Materials and Methods 

Data Acquisition 

We began assembling the wiring diagram by consolidating existing data from both published and 
unpublished sources. Using J. G. White et a/.'s The Mind of a Worm (MOW) Q as the starting point, 
we extracted wiring data from diagrams, figures, tables, and text (for example, see d Appendix A, 
pp. 118-122] on neuron AVAL/R). Connectivity of each neuron, its synaptic partner, and synaptic type 
(chemical, gap junction, neuromuscular) was manually entered into an electronic database. In the ventral 
cord of the worm, this level of synaptic specification was complicated by the fact that connections were 
recorded by neuron class. For example, bilateral neurons PVCL and PVCR were simply listed as PVC. 
We were able to assign proper connections to the appropriate left/right neuron by referring to White and 
coworker's original laboratory notebooks and original electron micrographs. In some cases, the number 
of synapses for a given neuron class in MOW differed from the sum of connections for the bilateral 
pairs in the notebooks and/or electron micrographs. The synaptic value of these neurons was determined 
by appropriating the value in MOW according to proportionality from the notebooks and/or electron 
micrographs. 

From here, we incorporated R. M. Durbin's data, which was applicable to the anterior portion of the 
worm, reconstructed from the animal N2U [8]. For neurons that projected beyond the nerve ring, only 
the anterior connections needed to be updated. Since data from MOW did not specify the location of 
synapses, integration proved difficult. For these neurons, we obtained positional information by cross- 
referencing Durbin's data against original electron micrographs and his handwritten annotations in White's 
laboratory notebooks. Only synapses located in regions addressed by Durbin were included. Connections 
in the middle and tail regions of the worm were mostly unaffected by these updates. 

With the advent of green fluorescent protein (GFP) reporters, researchers are able to visualize the 
neuroanatomy of individual neurons. Studies based on this technology mostly confirmed the electron 
micrograph reconstructions described in MOW. A few differences between GFP-stained neurons and 
White's work were observed [28]. Notably, the anterior processes of DVB and PVT could have been 
mistakenly switched in MOW [7]. Based on these findings, we reversed the connections for neurons DVB 
and PVT anterior to the vulva. 

Most published works have focused in the neck and tail regions of C. elegans, where the majority of 
neuron cell bodies reside. Reconstructions of neurons in the mid-body of the worm, on the other hand, 
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are scant and incomplete. From a combination of published works 0, flU, ifTOl . 1811 . we found that 
wiring data for 64 neurons had large gaps or were missing entirely. Sixty-one of these are motor neurons 
in the ventral cord. Two are excretory neurons (CANL/R) that do not appear to make any synapses. The 
remaining neuron, RID, is the only process in the dorsal cord that extends over the length of the animal. 

At the C. elegans archive (Albert Einstein College of Medicine), we uncovered a large number of 
reconstruction records in White et al.'s laboratory notebooks. These notebooks identified neurons by 
different color code labels depending on the animal, the location of the neurite (ventral or dorsal), 
and magnification of the electron micrograph. To ascertain the identity of the neurons, we relied on a 
combination of color code tables and comparisons of common anatomical structures between electron 
micrograph prints. In the end, we identified notes for full reconstructions of 24 of the aforementioned 
neurons. Partial connectivity data for the remaining 38 were also available where 22 neurons have 
partial/missing dorsal side connections and 6 neurons have partial ventral side connections. We checked 
the connections of all (both published and unpublished) neurons in the ventral cord against electron 
micrographs used by White and coworkers. Over 600 updates were made to the original notes and 
published reconstructions. Many of these updates were additions of previously missed neuromuscular 
junctions between ventral cord motor neurons and body wall muscles. 

While conducting this work, we found that a large section of the worm on the dorsal side, from 
just anterior to the vulva to the pre-anal ganglion, was never imaged at high power magnification with 
an electron microscope. This lack of electron micrographs was the reason why so many neurons were 
missing dorsal side reconstructions. Using original thin sections for the N2U worm prepared by White et 
al. , we produced new high power electron micrographs of this dorsal region. Due to the condition of the 
sections, only one of every 2-3 sections was imaged. These new electron micrographs extended nearly 
9/xm on the dorsal side. New dorsal side data for 3 neurons (DA5, DB4, DD3) were obtained from these 
electron micrographs. Resource constraints prevented us from covering the entire dorsal gap. 

From our compilation of wiring data, including new reconstructions of ventral cord motor neurons, we 
applied self-consistency criteria to isolate neurons with mismatched reciprocal records. The discrepancies 
were reconciled by checking against electron micrographs and the laboratory notebooks of White et 
al. Connections in the posterior region of the animal were also cross-referenced with reconstructions 
published in [10]. Reconciliation involved 561 synapses for 108 neurons (49% chemical "sends," 31% 
chemical "receives," and 20% electrical junctions). 



Giant Component for Random Networks 

For a random network with N neurons and probability p of a connection being present, if the constant 
c = Np > 1, then the size of the giant component is asymptotically normal with mean Na(c) and 
variance N(3{c) |[82l p. 120]. These quantities are given by 

7 7(1--) 

a ( c ) = 1 - I and /3(c) = ^, (14) 

c c(l — j) A 

where 

7 = "W(~). (15) 

and W(-) is the Lambert VF-function. If we take N to be 279 and p to be 514/ ( 27 2 9 ), then c = 3.698. Using 
the asymptotic approximation, the size of the giant component is distributed approximately normally with 
mean 271 and variance 9.22. Thus the probability of having a giant component of size 248, which is over 
7 standard deviations from the mean, is about 10 -14 . If a precise evaluation of this infinitesimal value is 
desired, large deviations techniques, rather than the asymptotic approximation may be more valid [83]. 
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To apply this method to the weakly connected component of a directed network, we are interested in 
the undirected network formed by adding a connection between two neurons if there is a connection in 
either direction. For a random directed network with probability q of presence of a directed connection, 
the probability of a connection existing in either direction is p = q 2 + 2q(l — q). Taking q to be 
2194/279/278 = 0.0283, p is 0.0558. Then for an undirected random network with N = 279 and 
the specified p, c is 15.56. Then the size of the giant component is distributed approximately normally 
with mean 279 and variance 0.0000487. Thus the probability of the giant weakly connected component 
containing all the neurons in such a random network is overwhelming. Again, large deviations techniques 
should be used to get a precise evaluation of the probability of being on the order of 10000 standard 
deviations away from the mean. 

Giant Component for Random Networks with Given Degree Distribution 

Consider the ensemble of random networks with a given degree distribution |[84l . For the gap junction 
network, the generating function corresponding to the measured degree distribution is 

G (x) = (279) [26 + 39x + 59x 2 + 43x 3 + 46x 4 + 23x 5 + 15x 6 + 5x 7 + 8x 8 + 4x 9 + 3X 11 + 2x 14 

+2 x 15 + x 24 + x 29 + x 34 + x 40 ] , 

with derivative 

G (x) = (275) [39 + 118x + 129x 2 + 184x 3 + 115x 4 + 90x 5 + 35x 6 + 64x 7 + 36x 8 + 33x 10 
+28x 13 + 30x 14 + 24x 23 + 29x 28 + 34x 33 + 40x 39 ] . 

Therefore G" (l) = 1 2 °^. The generating function G\ is then 

d(x) = (j^g) [39 + 118x + 129x 2 + 184x 3 + 115x 4 + 90x 5 + 35x 6 + 64x 7 + 36x 8 + 33x 10 + 28x 13 
+30x 14 + 24x 23 + 29x 28 + 34x 33 + 40x 39 ] . 

As shown in [84], the expected fraction of the network taken up by the giant component, S, is S = 
1 — Go(u), where u is the smallest non-negative solution to u = Gi(u). In our case, we find u = 0.043, 
and so S = 0.90. That is to say, one would expect the giant component to consist of 251 neurons. 

Using the computed S and G' (l), we can find the average component size excluding the giant 
component, which turns out to be 1.58. 

For the symmetrized chemical network, the generating function corresponding to the measured degree 
distribution is 

H (x) = (279) [2x + 6x 2 + 8x 3 + 6x 4 + 14x 5 + 14x 6 + 19x 7 + 20x 8 + 19x 9 + 20x 10 + 17X 11 + 18x 12 
+14x 13 + 9x 14 + 10x 15 + 9x 16 + 4x 17 + 9x 18 + 7x 19 + 3x 20 + 9x 21 + 8x 22 + 3x 23 + 4x 24 
+3x 25 + 2x 26 + 3x 27 + 2x 29 + x 31 + x 32 + 2x 33 + x 34 + x 36 + x 42 + x 48 + x 49 + 2x 50 
+x 51 + x 52 + x 53 + x 56 + x 83 + x 85 ] , 

with derivative 

H' (x) = (279) [2 + 12x + 24x 2 + 24x 3 + 70x 4 + 84x 5 + 133x 6 + 160x 7 + 171x 8 + 200x 9 + 187x 10 
+216X 11 + 182x 12 + 126x 13 + 150x 14 + 151x 15 + 68x 16 + 162x 17 + 133x 18 + 60x 19 + 189x 20 
-f 176x 21 + 69x 22 + 96x 23 + 75x 24 + 52x 25 + 81x 26 + 58x 28 + 31x 30 + 32x 31 + 66x 32 
+34x 33 + 36x 35 + 42x 41 + 48x 47 + 49x 48 + 100x 49 + 51x 50 + 52x 51 + 53x 52 + 56x 55 + 83x 82 
+85x 84 l . 
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Therefore H' (l) = ^§-- The generating function H\ is then 

Hi(x) = (3^9) [2 + 12s + 24x 2 + 24x 3 + 70x 4 + 84x 5 + 133x 6 + 160x 7 + 171x 8 + 200x 9 + 187x 10 
+216X 11 + 182x 12 + 126x 13 + 150x 14 + 151x 15 + 68x 16 + 162x 17 + 133x 18 + 60x 19 + 189x 20 
+176x 21 + 69x 22 + 96x 23 + 75x 24 + 52x 25 + 81x 26 + 58x 28 + 31x 30 + 32x 31 + 66x 32 
+34x 33 + 36x 35 + 42x 41 + 48x 47 + 49x 48 + 100x 49 + 51x 50 + 52x 51 + 53x 52 + 56x 55 + 83x 82 
+85x 84 ] . 

The expected fraction of the network taken up by the giant component, S, is S = 1 — Hq(u), where u 
is the smallest non-negative solution to u = H\{u). Here u is found to be 0.00051, and so S = 0.999996. 
That is to say, one would expect the giant component to consist of 278.9990 neurons. 



Path Length for Random Networks with Given Degree Distribution 

Continuing from the previous subsection, we find the derivative of the generating function G± for the 
gap junction network to be 

G[(x) = (-y [118 + 258x + 552x 2 + 460x 3 + 450x 4 + 210x 5 + 448x 6 + 288x 7 + 330x 9 + 364x 12 
+420x 13 + 552x 22 + 8 1 2x 27 + 1122x 32 + 1560x 38 ] . 

Thus G[(l) = Letting z t = G' (l) = and z 2 = = ^f, it is shown in [84, (53)], 

that the expected path length is 

ln[(N-l)(z 2 - Zl ) + zf] -Inzf 
L = — L - — - — — J - = 3.05. (16) 

In [z 2 /zi\ 



Fitting Tails of Distributions 

To find functional forms of the tails of various distributions, we follow the procedure outlined in 
Il44l . For the candidate functional forms — say, the power law p(d) ~ dr 1 and the exponential decay 
p(d) ~ exp(— Xd) — we perform the following steps. First, we find the optimal parameter of the fit by 
maximizing the log-likelihood and the optimal starting point of the fit by minimizing the Kolmogorov- 
Smirnov statistic. Second, we evaluate the goodness of fit by calculating the p-value that the observed 
data was generated by the optimized distribution using p > 0.1 as a criterion for plausibility. Finally, if 
several distributions pass the p-value test we compare their log-likelihoods to find the most probable one. 



Circuits in Eigenmodes 

Let us bound the probability of finding an eigenmode that comprises a random set of neurons. Let N 
be the number of neurons in the network being analyzed. Let Ki be the number of neurons that appear 
strongly in the ith eigenmode and let K = maxj Ki. Furthermore let M be the number of neurons in 
the random set, which one might endeavor to investigate as a putative functional circuit derived from an 
eigenmode. 

Now go through each eigenmode and add to a list all possible unordered M-tuples of strong neurons. 
Even if all of these are unique, the size of the list is upper-bounded by ^i=i (m) which itself is 
upper-bounded by N (^) . 

Additionally, we can compute the number of all unordered M-tuples of neurons. This number is (^) . 
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Thus, if a random set of neurons was selected from all possible sets of neurons, the probability p that 
there would be an eigenmode containing all of them is upper-bounded as 

< Elli (m) < HSl = NK\ M\(N-M)\ _ jf(if-l)...(jf-M + l) 
P ~ (jj) ~ (jj) (2T-M)!M! iV! (JV-l)(iV-2)---(iV-M+l) - iV M - r 

Suppose we are interested in putative functional circuits of size M = 6 in the giant component of 
the gap junction network, which has N = 248 and from Figure |S2l take K = 20. Then even the loosest 
upper-bound yields 

K M 2D 6 

P < ~ior^ = f = 6.8 x 10~ 5 , 

* - N m-i 248 5 

and so finding a random set of neurons in an eigenmode is unlikely. 

Suppose we know L functional circuits of size M through molecular biology and want to know the 
probability of at least one of them appearing in the eigenmodes by chance. By the union bound (Boole's 
inequality), this probability is less than pL. If we take L = 20 and M = 6, the probability of a known 
functional circuit appearing in the eigenmodes by chance is less than 1.4 x 10 -3 for the giant component 
of the gap junction network. 

Gap Junction— Chemical Synapse Likelihoods 

The likelihood ratios shown in Figure [9] are the following quantities, empirically estimated from either 
all neuron pairs or pairs with a GABAergic presynaptic neuron. The first is 

Prfchem. absent|gap absent] 



The second is 



and the third is 



Pr[chem. absentjgap present]' 

Prfchem. unidirectional | gap absent] 
Prfchem. unidirectional | gap present]' 

Pr[chem. bidirectional | gap absent] 
Pr[chem. bidirectional | gap present] 



Degree Correlation Coefficients 

Table U shows the correlation coefficients between neuron degree sequences. The average correlation 
coefficients of randomly permuted degree sequences from 10000 trials are also shown for comparison. The 
standard deviation is also shown since the distributions of the three randomized correlation coefficients 
were all nearly symmetric about zero. 



e-Pseudospectrum Computation 

We used the MATLAB package EigTool IIB31 to compute pseudospectra. 

MATLAB Code and Data 

Note that MATLAB code for computing several network properties is available at 
|http : //mit . edu/lrv/www/elegans/| This collection of software may be used not only to re- 
produce most of the figures in this paper, but also for future connectomics analyses. 

The collected data is available from the WormAtlas 11221 as well as from the same website as the 
MATLAB code. 
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Appendix A 
Algorithm for Directed Network Drawing 

To visualize a directed neuronal network we modify an approach suggested in 11361 , IT381 . In this 
approach, the vertical and the horizontal coordinates are chosen independently. The arrangement of 
neurons along the vertical axis conveys information about the directionality of the signal flow in the 
network and the arrangement of neurons along the horizontal axis or axes conveys information about the 
strength of connectivity regardless of directionality. 

To find the vertical coordinate, z, we try to arrange the neurons so that for every synaptically connected 
pair of neurons, the difference in z between a presynaptic neuron i and a postsynaptic neuron j is as 
close to one as possible. Specifically, we minimize the following energy function: 

1 n 

E =2 22 W *i & ~ z i ~ s § n ^ " A ^ (17) 

of the connectivity matrix Aij, which is the sum of the gap junction and chemical connectivity matrices, 
and the symmetrized connectivity matrix Wij, which satisfies Wy = (Aij + Aji)/2. By setting the 
derivative of this expression to zero, we find: 

Lz = b, (18) 

where 6, = Y^=i WySgn(^4y — A^) and the Laplacian L = D — W is defined in terms of a diagonal 
matrix D that contains the number of synaptic terminals on corresponding neurons, 

n 

Dij = SijJ^Wik- (19) 

k=l 

A unique solution to this equation can be found by using the pseudoinverse. 

To find the horizontal coordinates, we use the Laplacian, L, normalized by the number-of-terminals 
matrix D, 

Q = D- 1/2 LD~ 1/2 . (20) 

The eigenmodes corresponding to the second and third lowest eigenvalues of Q are denoted v 2 and v$. 
Then, the horizontal coordinates are 

x = D~ 1/2 v 2 and y = D~ 1/2 v 3 . (21) 

This method produces an aesthetically appealing drawing because each neuron is placed in the weighted 
centroid of its neighbors. Thus strongly coupled neurons tend to be colocated. 
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Appendix B 
Algebraic Form of Survival Functions 



Here we consider several commonly encountered distributions and their survival functions. If a dis- 
tribution were to follow a power law, p(d) ~ d~ a , then the survival function (under a continuous 
approximation) also follows a power law: 

oo 

P ( d ) = Y^k- a -d-^-V. (22) 

k=d 

Similarly, if a distribution follows an exponential decay, p(d) ~ e~ d l K , then the survival function also 
has an exponential decay, with the same exponent: 

oo 

P(d) = e~ k/K ~ e~ d / K . (23) 

k=d 

If a distribution were to follow the (continuous) stretched exponential distribution, p(d) ~ (d/ /3) 7_1 e _ ( d /^ 7 , 
then the survival function would have a decay given by a stretched exponential function with the same 
stretch factor 7: 

00 

P(d) = ^(k/py-^e-WW ~ e-WW. (24) 
k=d 

Appendix C 

Further Spectral Properties of the Gap Junction Network 

The spectral radius of the Laplacian plays a significant role in the performance of linear systems with 
dynamics that are slightly different from charge equilibration, but which have been used to describe the 
synchronization of networks of oscillators and the operation of distributed control systems in engineering 
[Ml, ll87ll . The spectral radius is the largest eigenvalue and is denoted by Aat for networks of size N. 
For these dynamics, the ratio of the spectral radius and the algebraic connectivity, Aat/A2 determines the 



rate of convergence of synchronization. From Figure |4(a)[ it may be computed that the eigenratio for the 
giant component of the actual gap junction network is 1026. 
There is a general lower bound for the eigenratio 



A 2 «min 

where d mSLX is the maximum degree of all the neurons in the network. There are networks which achieve 
this bound. For a network with maximum degree 40 and minimum degree 1, as in the giant component 
of the gap junction network, this bound is 41. We see that the eigenratio for the actual gap junction 
network is 1026, off from the optimal. 

Another quantity that often arises in discussions of signal propagation in networks is the magnification 
coefficient, c [ 89 ]— [92 ]. Networks that have large magnification coefficients transmit signals quickly. The 
magnification coefficient is difficult to compute, but can be approximated by the algebraic connectivity 
through an unexpected connection between local connectivity properties and spectral properties. A large 
algebraic connectivity implies a large magnification coefficient. In particular, there are upper and lower 
bounds that relate the two lf89l . 

c > + 2A 2 , (26) 

^max 
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and 



c 2 

A 2 > — + 2c 2 . (27) 
The algebraic connectivity may be compared to a general upper bound [93]: 

A 2 < , (28) 

where d m j n is the minimum degree of all the neurons in the network. There are classes of networks, 
called Ramanujan graphs that have constant degree and that can get close to the bound. For a network 
with 248 neurons, 511 connections, and minimum degree 1, as in the giant component of the gap junction 
network, the bound is 1.00. For the gap junction giant component, the algebraic connectivity is 0.12, so 
we see that it is not very close to the upper bound, but is not too far either; the algebraic connectivity 
is reasonably large. The main cause for deviating from the bound is non-constant degree distribution. 

Since the algebraic connectivity is fairly large, the gap junction network also has a fairly large 
magnification coefficient. 

Appendix D 

ElGENDECOMPOSITION 

Physical systems are often represented by linear, constant-coefficient differential equations. Differential 
equations provide an implicit specification of the system, giving the relationship between input and output, 
rather than an explicit expression for the system output as a function of the input. After specifying initial 
conditions, differential equations can be solved to find explicit expressions for the output. 

Dynamical systems that can store energy in only one form and location are called first-order, since the 
equation describing time evolution can be written only in terms of a single variable and its first derivative. 
Storing energy is a form of short-term memory. For a single state variable Vi, a canonical first-order, 
linear, constant-coefficient differential equation is 

at 

where r and La are fixed constants and M(t) is some signal. 

The natural (unforced) response of a system corresponds to M{t) = and is completely determined 
by the system's eigenvalue. In particular, solving 

with initial condition Vi(t = 0) = Vq, yields 

v t (t) = y e ( - Wr) ', 

where —Lu/t is the eigenvalue. 

The forced response occurs when some exogenous perturbation is applied to the system. For example 
if a scaled step function Mou(t) is applied, then the differential equation 

T ^. + L u V i {t)=M u(t) 
at 

with initial condition Vi(t = 0) = Vq has solution 

M 



-Liii 



V - L 



11 



e (-Z„/r)t ^ t>Q 
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The response of a first-order system to a unit impulse is identical to its natural response; the impulse 
generates the initial condition in such a short time that it has no other effect on the system. That is, the 
system is jarred to the initial position by the impulse. 

Generally when a forcing function is applied to a linear constant-coefficient dynamic system, the 
response will consist of the superposition of the forced response (a modification of the input signal) and 
the natural response governed by the system's eigenproperties. 

Thus far, we have considered a single state variable Vi(t), but in neuronal networks we actually 
have a vector of states, V(t) = [Vx(t) V 2 {t) • • • Vjv(t)] > governed by a system of linear constant- 
coefficient differential equations. A canonical form is 

J,] '-L 11 V 1 (t) + L 12 V 2 (t) + --- + L 1N V N (t) = M 1 (t) 



dt 
_dV 2 {t) 

dt 



+ L 2l Vx{t) + L 22 V 2 (t) + ■■■ + L 2N V N (t) = M 2 (t) 



r dVN(t)_ + LmVi ^ + Lm v 2 {t) + ■■■ + L NN V N {t) = M N {t) 



dt 

which can be written in matrix-vector form as 

dVit) 

T 

dt 



+ LV(t) = M(t). 



The natural response of such a system with initial condition V(t = 0) = Vq is the vector 

V(t) = Voe { ~ L/T)t . 

Although this is in principle the solution to the system of differential equations, it is difficult to examine. 
Study of system behavior is complicated by the fact that each of the equations is coupled to the others 
through the off-diagonal elements of L. It would be desirable to find a new coordinate system in which 
all equations are decoupled (such that the coefficient matrix is diagonal). 
A vector v is called an eigenmode of a matrix L if it satisfies 

Lv = Xv 

for some number A, which is called the eigenvalue. Decomposing the coefficient matrix into its eigen- 
decomposition0 



L = [vi v 2 



VN 



Ai 






A 2 










Ajv 



[vi v 2 



VN 



-1 



allows us to write the natural response as 

V(t) 



N 

E 

i=l 



(-Vr)t 



Note that not all matrices have an eigendecomposition. Instead, the Jordan decomposition should be used for these non- 
diagonalizable matrices |94|. The three matrices we consider, L, A T , and $ are diagonalizable and so the eigendecomposition 
is identical to the Jordan decomposition. 

Another decomposition that has been proposed for use in systems neuroscience is the Schur decomposition |95|. Since the 
gap junction network is undirected, the Schur decomposition is also identical to the eigendecomposition. For the chemical and 
combined networks, the Schur modes may provide additional insights, but we do not consider them in this work. 
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where ctj is the projection of the initial condition vector Vq onto Uj. 

The essential idea of the eigenmode decomposition is that the natural response of the system can 
be viewed as the superposition of a number of distinct types of dynamics — the eigenmodes — each one 
associated with a particular natural frequency of the system. The natural frequencies, — \/t, of the 
system are determined by the eigenvalues Aj of L. Each mode involves excitation of one and only one 
natural frequency of the system. 

If an eigenmode is real, then the dynamics associated with the solution can be described by a straight 
line in state space. The system moves in the direction of the eigenmode. For example, moving in the 
direction of the eigenmode [+1 —1 • • • 0] would equalize the values of V% and V2 but 
not affect V3, . . . , Vjy. A more complicated eigenmode would involve all state variables that are non-zero. 

Beyond their simple geometric interpretation in state space, the eigenmodes also have a simple repre- 
sentation as time functions, since each one involves a single exponential rather than a mixture of several 
exponentials with different exponents. The exponent — \/t determines how quickly the system response 
in the direction of eigenmode Vi decays. For fixed r, the larger the eigenvalue Aj, the more quickly the 
eigenmode decays. 

The forced response of a network proceeds in the same way as the forced response of a scalar system. 
Further details on linear system analysis with eigenmodes can be found, e.g. in the textbooks l96ll . |97l . 
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Supplemental Material 



table si 

Connected components of the gap junction network. Note the single giant component and the large 

number of disconnected/isolated neurons. 



Giant Component (248 neurons) 


ADAL/R 


ALNL 


AVG 


DD01-05 


PDA 


PVR 


RIVL/R 


SABVL/R 


URYVL/R 


ADEL/R 


AQR 


AVHL/R 


DVA 


PDB 


PVT 


RMDDL/R 


SDQL/R 


VA01-12 


ADFL/R 


AS01-11 


AVJL/R 


DVB 


PDEL/R 


PVWL/R 


RMDL/R 


SIADL/R 


VB01-11 


ADLL/R 


ASGL/R 


AVKL/R 


DVC 


PHAL/R 


RIBL/R 


RMDVL/R 


SIAVL/R 


VC01-05 


AFDL/R 


ASHL/R 


AVL 


FLPL/R 


PHBL/R 


RICL/R 


RMED 


SIBDL/R 


VDOl-10,13 


AIAL/R 


ASIL/R 


AVM 


IL1DL/R 


PHVL/R 


RID 


RMEL/R 


SIBVL/R 




AIBL/R 


ASKL/R 


AWAL/R 


IL1L/R 


PLML/R 


RIFL/R 


RMEV 


SMBDL/R 




AIML 


AUAL/R 


AWBL/R 


IL1VL/R 


PQR 


RIGL/R 


RMFL 


SMBVL/R 




AINL/R 


AVAL/R 


BAGL/R 


IL2L/R 


PVCL/R 


RIH 


RMGL/R 


SMDDL/R 




AIYL/R 


AVBL/R 


CEPDL/R 


LUAL/R 


PVM 


RIML/R 


RMHL/R 


SMDVL/R 




AIZL/R 


AVDL/R 


CEPVL/R 


OLLL/R 


PVNL 


RIPL/R 


SAADL/R 


URBL/R 




ALA 


AVEL/R 


DA01-09 


OLQDL/R 


PVPL/R 


RIR 


SAAVL/R 


URXL/R 




ALML/R 


AVFL/R 


DB01-07 


OLQVL/R 


PVQL/R 


RIS 


SABD 


URYDL/R 





First Small Component (2 neurons) 

ASJL/R 



Second Small Component (3 neurons) 

HSNL/R PVNR 

Neurons with no gap junctions (26 neurons) 

AIMR ASEL/R BDUL/R IL2DL/R PLNL/R RIAL/R URADL/R VD11-12 
ALNR AWCL/R DD06 IL2VL/R PVDL/R RMFR URAVL/R 



TABLE S2 

(A) Number of gap junction contacts between different neuron categories. (B) Percent of gap junctions 

ON NEURONS OF THE ROW CATEGORY THAT CONNECT TO NEURONS OF THE COLUMN CATEGORY. 



A 


Sensory 


Inter- 


Motor 


Sensory 


108 


119 


26 


Inter- 


119 


368 


342 


Motor 


26 


342 


324 



B 


Sensory 


Inter- 


Motor 


Sensory 


42.7% 


47.0% 


10.3% 


Inter- 


14.4% 


44.4% 


41.3% 


Motor 


3.8% 


49.4% 


46.8% 
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TABLE S3 

Comparison of clustering coefficient and characteristic path length of the giant component of the C. 

elegans GAP JUNCTION NETWORK AND SEVERAL OTHER NETWORKS THAT HAVE BEEN CLASSIFIED AS SMALL WORLD 
NETWORKS. THE CLUSTERING COEFFICIENT OF AN EQUIVALENT ERDOS-RENYI RANDOM NETWORK IS INDICATED IN 
PARENTHESES. THIS IS CALCULATED USING THE WATTS AND STROGATZ APPROXIMATIONS TO L AND C BY FINDING 

a«iexp(Mp). 



Network 


N 


C (CV) 


L 


Giant component of gap junction network 


248 


0.21 (0.014) 


4.52 


Analog electronic circuit |52| 


329 


0.34 (0.019) 


3.17 


Class dependency graph of Java computer language |98| 


1376 


0.06 (0.002) 


6.39 


Film Actors |14] 


225226 


0.79 (0.00013) 


3.65 


Power Grid |14| 


4941 


0.080 (0.00032) 


18.7 



TABLE S4 

(A) Number of chemical synapse contacts from row category to column category. (B) Percent of 

SYNAPSES IN ROW CATEGORY THAT SYNAPSE TO COLUMN CATEGORY. 



A 


Sensory 


Inter- 


Motor 


Sensory 


474 


1434 


353 


Inter- 


208 


1359 


929 


Motor 


30 


275 


1332 



B 


Sensory 


Inter- 


Motor 


Sensory 


21.0% 


63.4% 


15.6% 


Inter- 


8.3% 


54.5% 


37.2% 


Motor 


1.8% 


16.8% 


81.4% 



TABLE S5 

Strongly connected components of the chemical network. Note the single giant component and the 

large number of isolated neurons. 



Giant Component (237 neurons) 


ADAL/R 


ALNL/R 


AVFL/R 


CEPVL/R 


LUAL/R 


PVM 


RIH 


RMHL/R 


URYDL/R 


ADEL/R 


AQR 


AVG 


DA0 1-06,09 


OLLL/R 


PVNL/R 


RIML/R 


SAADL/R 


URYVL/R 


ADFL/R 


ASOl-06,09,11 


AVHL/R 


DB0 1-04,07 


OLQDL/R 


PVPL/R 


RIPL/R 


SAAVL/R 


VA0 1-06,08-09, 11-12 


ADLL/R 


ASEL/R 


AVJL/R 


DD0 1-02,05 


OLQVL/R 


PVQL/R 


RIR 


SABD 


VB0 1-06,08- 11 


AFDL/R 


ASGL/R 


AVKL/R 


DVA 


PDA/B 


PVR 


RIS 


SDQL 


VC01-05 


AIAL/R 


ASHL/R 


AVL 


DVC 


PDEL/R 


PVT 


RIVL/R 


SMBDL/R 


VDOl-03,05-06,08, 10-13 


AIBL/R 


ASJL/R 


AVM 


FLPL/R 


PHAL/R 


PVWL/R 


RMDDR 


SMBVL/R 




AIML/R 


ASKL/R 


AWAL/R 


HSNL/R 


PHBL/R 


RIAL/R 


RMDL/R 


SMDDL/R 




AINR 


AUAL/R 


AWBL/R 


IL1DL/R 


PLMR 


RIBL/R 


RMDVL 


SMDVL/R 




AIYL/R 


AVAL/R 


AWCL/R 


IL1L/R 


PLNL 


RICL/R 


RMED 


URADL/R 




AIZL/R 


AVBL/R 


BAGL/R 


IL1VL/R 


PQR 


RID 


RMEV 


URAVL/R 




ALA 


AVDL/R 


BDUL/R 


IL2L/R 


PVCL/R 


RIFL/R 


RMFL/R 


URBL/R 




ALML/R 


AVEL/R 


CEPDL/R 


IL2VL/R 


PVDL 


RIGL/R 


RMGL/R 


URXL/R 





Small Component (2 neurons) 

RMDVR RMDDL 



Isolated neurons in chemical network (40 neurons) 

AIM! DA07-08 DVB PIJVIL RMEL/R SDQR SIAVL/R SIBVL/R VB07 

ASIL/R DB05-06 IL2DL/R PLNR SABVL/R SIADL/R SIBDL/R VA07,10 VD04,07,09 

AS07,08,10 DD03-04,06 PHCL/R PVDR 
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Fig. SI. Geodesic distance distributions. |(a)| Giant component of gap junction network. |(b)| Giant component of chemical 
network. |(c)[ Giant component of combined network. 



TABLE S6 

SOME STRUCTURAL PROPERTIES OF THE C. elegans GAP JUNCTION NETWORK, RANDOMLY EDITED NETWORKS (E gap ), AND 

THE AY NETWORK lfl3ll . 





C. elegans 


AY's C. elegans |13| 


Egap 






454 


177 ±18.5 


giant component size 


248 


253 


261 ±3.41 


giant component pathlength 


4.52 


4.71 


4.09 ± 0.078 


giant component clust. coef. 


0.21 


0.23 


0.14 ±0.011 



TABLE S7 

Some structural properties of the C. elegans chemical network, randomly edited networks (E C hem), and 

THE AY NETWORK lfl3ll . 





C. elegans 


AY's C. elegans |13| 


Echcm 






3546 


638 ± 33.2 


weak giant component size 


279 


279 


279 ± 0.07 


strong giant component size 


237 


239 


267 ±3.19 


strong giant component pathlength 


3.48 


3.99 


3.12 ±0.028 


strong giant component clust. coef. 


0.22 


0.20 


0.16 ±0.006 
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Fig. S2. Eigenmodes of Laplacian for giant component of gap junction network. 
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Fig. S3. Fastest eigenmodes of Laplacian for giant component of gap junction network. Eigenmodes corresponding to 
A248, A247, • • • , A234 are shown. The eigenmodes are labeled with neurons that take value above a fixed absolute value threshold. 
Neurons with negative values are in red, whereas neurons with positive values are in black. 
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